function randu(iseed)
!-----------------------------------------------------------------------
!c     generates randu pseudo uniform [0,1] random number generator.
!c     requires a seed (iseed) which is updated and returned.
!
! so, at each exit from this subroutine seed is changed. To simulate a
! sequence, start seed at a number, repeatedly call this subroutine
! to generate vars. Then, when simulating again, make sure seed is
! reset to the initial value so you get the same sequence.
!
!c     [lewis,goodman & miller,ibm systems journal,v.8#2 1969,pp.136-145]
! this paper is saved in S:\Hajivass_fort_code
!
! don't change iseed between calls
! #'s for iseed should be big: 12345678, 855998726,745681489,
! 506104362,236686234,1912615462,864268549,13034519
!c-----------------------------------------------------------------------
INTEGER, intent(inout)   ::iseed
DOUBLE PRECISION   ::seed,xb,xm,randu
INTEGER          :: xseed

!data xb,xm/16807.d0,2147483647.d0/
xb=16807.0D0
xm=2147483647.0D0
seed=dble(iseed)
xseed=mod(xb*seed,xm)
randu=xseed/xm
iseed=xseed
end function randu

